Radiation spectra of laser-driven quantum 
relativistic electrons 



Guido R. Mocken and Christoph H. Keitel 

Max- Planck- Institut fur Kernphysik, Saupfercheckweg 1, D-69117 Heidelberg, 

Germany 

Theoretische Quantendynamik, Physikalisches Institut, Universitdt Freiburg, 
Hermann-Herder-StrajJe 3, D-79104 Freiburg, Germany 



Abstract 

A procedure to calculate the radiation spectrum emitted by an arbitrarily prepared 
Dirac wave packet is developed. It is based on the Dirac charge current and classical 
electrodynamic theory. Apart from giving absolute intensity values, it is exact in 
terms of relativistic retardation effects and angular dependence. We employ a laser 
driven free electron to demonstrate the advantages of our method as compared to 
traditional ones that merely rely on the Fourier transform of the dipole operator's 
expectation value. Classical reference calculations confirm the results obtained for 
the low-frequency part of the spectrum, especially in terms of the observed red-shifts, 
which clearly deviate from non-relativistic calculations. In the high-frequency part 
of the spectrum, we note appreciable deviations to the purely classical calculations 
which may be linked to quantum averaging effects. 
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1 Introduction 



The appearance of harmonics in the light scattered at free electrons is known 
both in theory [l[ 0- 0, 0] as well as in experiments 0, S for a rather long 
time now. Atomic systems 0,0, Ell 11, HI 3] i n highly intense laser fields have 



turned out to produce especially high orders, the so-called High Harmonics 
jlil EH 3]- In this case, the experimentally observed spectra [13,13,11^ are 
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usually explained in terms of the so-called recollision model [20]. Therefore, 
in the context of quantum mechanical simulations such as those shown in 
apart from the obvious observables such as the expectation value 



of the position or momentum operator and the spatial probability density 
distribution, there is a lot of interest in calculating the radiation spectrum as 
well. 



The numerical methods that were described in our previous work 22] pro- 
vide the Dirac wave function as a function of time and position, where both 
variables are discretized and limited to certain regions. There exists a well- 
established procedur e |24f that derives the emission spectrum from such a 
given data set (e.g. [12J, |25(). However, this one was developed for the non- 
relativistic Schrodinger theory, and contains further approximations. The ques- 
tion arises whether, in the highly relativistic regime, where non-relativistic 
quantum theory breaks down, it is still advisable to employ a partly non- 
relativistic method for the calculation of the quantum emission spectrum. 



In this paper, we present a procedure that takes into account all relativistic 
spectral effects. It is based on the Dirac charge current and the classical theory 
of electromagnetism. Apart from giving absolute intensity values, it is exact 
in terms of relativistic retardation effects and angular dependence. To take 
advantage of the new scheme, one needs to obtain the time-resolved Dirac 
wave function of some physically interesting scenario in the first place. We 
made use of our advanced two-dimensional Dirac split-operator code 22|| and 
considered the motion of a free electron in a laser field to put the new 
scheme to the test. Classical calculations, which should not deviate much from 
the quantum calculations for this particular case as far as the lower part of 
the spectrum is concerned, serve as a reference and prove the correctness of 
the new scheme. However, significant deviations from the established method 



24j are revealed. In addition to that, the higher order harmonics are averaged 



out. The delocalized quantum charge distribution is likely to be responsible 
for this effect. 



This remainder of this paper is organized as follows: Section 2 briefly reviews 
the formulas required to evaluate the radiation spectrum of a classical point 
charge and their non-relativistic limit. Section 3 briefly discusses the estab- 
lished standard methods, then turns to the derivation of our new method. In 
the end, the non-relativistic limit of the latter is evaluated to show that it 
matches the former. This section also discusses the solution of some numerical 
problems that are involved. Section 4 discusses the results obtained with our 
method and relates them to those obtained with the standard methods and 
classical formulas. We draw our conclusions in section 5 and provide the clas- 
sical equations of motion in appendix A and the Dirac equation in appendix 
B, as well as some remarks and references concerning the way that we solved 
both problems numerically. 
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2 Classical radiation spectrum 



The physical quantity of interest here is the differential amount of energy 
(d 2 W) emitted into a certain interval of frequency (du) and solid angle (dft), 
i.e. the quantity -f^fj. 



2. 1 Spectrum of a charged classical point particle in relativistic motion 



From [27[ we obtain a formula for the far-field radiation of a classical charged 
point particle which is ready for application: 



d 2 W _ . q 2 



(X, UJ 



duidVt ' An 2 c 



+9° n x ( (n - (3) x (3 

dt — ^ L_^ e H^) 



where R(t) = \x - y(t)\ , n(t) = 0(f) = ® = I^(t). (1) 

Here, y(f) represents the particle's trajectory and x the observer's position. 
The particle's charge is labelled q, and c is the speed of light. The integral 
/ dt(. . .) is replaced by a sum At(. . .) which is carried out for all values u 
that are of interest by using the discrete trajectory data for y, y, y, that is 
available from the numerical integration of the classical equation of motion 
(A.5). 



2.2 Non-relativistic limit 



If one omits relativistic corrections (\/3\ 1) and accounts (with respect to 
the retardation R) only for particle displacements in the direction towards the 
observer that are very small (x ■ y ~ 0) as compared to the distance of the 
fixed observer \x\, then R is constant, 



R = \x — y\ =' x ' \x\ — x ■ y + O ( ] ~ const. — x ■ y Is const. (2) 




and instead of the integral in equation (1), the result is just the Fourier trans- 
form of the acceleration component a± that is perpendicular to the momentary 
direction of observation n: 
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d 2 W 
dudVL 



lx,u) 



4:7T 2 C 

2-kc 



dt n x [fix (3) e 



2vrc 



1 

2lr 



dt(3 ±{t) e iut 



-/5±(-w) 



2 g 2 



2nc 



(3) 



Using /3± D =' — one finally obtains: 



d 2 W 
dudQ 



2ttc 3 



a>±(u) 



(4) 



3 Quantum mechanical treatment 



The usual way (e.g. 12[ 25J) to calculate the emission spectrum is the follow- 
ing: Starting from the classical non-relativistic formula (4) the spectrum for 
a quantum mechanical system is evaluated using the corresponding quantity 
that is based on the acceleration expectation value (a): 

d 2 W q 2 2 

= WM ( 5 ) 



dujdVt ' 2irc 3 

For a more precise derivation and in addition to that a detailed discussion of 
the role of the "coherent" part (see above) and the "incoherent" part (not dis- 
cussed here) of the spectrum, see ji^, . The next step |24| is the application 
of the Ehrenfest theorem, i.e. the replacement 

/ — a Ehrenfest / Lorentz- / . Q -* ft\ //-,\ 

m ^ m = . F ) f = {QE+-VXB), (6) 

Newton force C 

which is easy to evaluate if the magnetic field B is neglected (the so-called 
dipole approximation). However, it is the magnetic field effects and the for- 
ward drift caused by them which have to be taken care of in this work with 
an emphasis on relativistic velocities. This precludes equation (4) from be- 
ing the basis of any highly-relativistic theory. In a weakly-relativistic work 
such as [30], the starting point remains equation (4), however, the relativis- 
tically correct equation (A. 3) is employed instead of equation (6). Utilising 
several essential approximations, the authors obtain an expression that can 
be evaluated numerically. Those approximations are no longer valid in the 
more strongly relativistic cases considered here, where even for classical calcu- 
lations, the substitution of equation (1) by (4) can already lead to noticeable 
frequency shifts in the spectrum. 

Unfortunately, in a quantum mechanical context, equation (4) is not easily 
replaced by the relativistically correct variant (1) and a calculation of the 
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required expectation values in the framework of the Dirac theory. On the 
one hand there are fundamental reasons, such as the difficulty to define a 
reasonable velocity operator in the Dirac theory 1 , on the other hand one 
needs to take into account certain practical considerations 2 . Therefore we 
have dropped all corresponding considerations in favour of the following far 
more promising idea. 



We start by calculating the emission spectrum of a spatially distributed clas- 
sical charge current, instead of the one of a point charge. The following deriva- 
tion differs slightly from those found in the literature [27] in that the one and 
only approximation made is the far-field expression for the electric field of a 
charge current. Otherwise it is exact. The transition to quantum mechanics 
is carried out at the very end of the calculation by substituting a suitable 
expression for the charge current. 



3. 1 Fully relativistic spectrum of a spatially distributed charge current 



If we call the magnitude of the Poynting-vector (energy per area and unit 
time) l^l, its direction ft, and the electric and magnetic field E and B, then 
we obtain 



S = ±(E x B) §= ^ xg ^-(Ex(nx E)) H ^ =0 ±E 2 n 
^\S\ = ^-E 2 (7) 

We assume the radiation detector to be located at position f = rf and to 
be activated during the time interval [to^i]- Let W(i) be a function 3 that 
describes the detection interval: 

m) = iKt)>o v*eMi], (8) 

otherwise. 

If we designate energy by W, power by P and angular frequency by u), then 
the energy per unit solid angle Q collected by the detector at position f in the 
time interval of interest [toj^i] is 



1 See the discussion of Zitterbewegung and the Foldy-Wouthuysen transform in 
0. 

2 Operator functions that contain quotients and square roots, and in which mixed 
position and momentum operator dependent terms occur (symmetrized as neces- 
sary), can be evaluated only with disproportionately large numerical effort. 

3 Later on, we will restrict ourselves to the simple case h(t) = 1. 
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Ut W(t)-(t) 

— oo 



dP=\S\r 2 dn 
(7) 



+00 



Jdt W(t)£r\E(t,f)) 2 - 

-00 

Here, we insert the Fourier transform of the electric field and its inverse: 



(9) 
(10) 



+00 

E(co, r) = -= Jdt E(t, f)e-^ ^ T (E {t , n ) 



(u>,r) 



-00 

Koo 



E(t, r) = -L jdu E(uj, ^ T' 1 (e^) ^ 



'ID 



(12) 



Since ^WE is real, we have T ( VWE) = T ( VWE) j and therefore one 



finally arrives at 

dW 

The integrand of (13) yields 

d 2 W 



'(-«) 



(13) 



dudQ, 



c 

27 



(14) 



The electric field itself is given by the time derivative of a certain vector poten- 
tial Af as obtained in the usual far field approximation for a classical charge 
current j(t',f') that is essentially localized at the origin of the coordinate 
system: 



^• r) sAwt A ' {t - n - where (15 » 

A f (t, r) D = ^- J d 3 r' r x {j(ct - r + r ■ r' , f') x f } + O Q^j . (16) 



Here, f is a constant unit vector that points towards the observer who is sit- 
uated at a fixed distance r, i.e. at position f— rr. If we insert (16) into (15), 
convert the result to the frequency domain using (11) and, while keeping in 
mind (8), insert everything into (14), we finally obtain, neglecting contribu- 
tions that vanish as O (jj for large r 
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d 2 W 
dudVl 



r =ct 



4tt 2 c 3 



cti 



dr ld 3 r' 



cto R 3 



r x 



_d_ 

dr 



J*(r — r + f 



f f') 



x r) 



(17) 



In the above equation and some more to follow the symbol "(g)" is just a 
multiplicative line continuation marker, which we have to distinguish from the 
cross product "x" . Note that so far, the calculation is an entirely classical one. 
We are aiming at replacing the classical current by the quantum mechanical 
one. Since our numerical Dirac code delivers ip and hence j at fixed, equally 
spaced times, whereas equation (17) requires very different times depending 
on position r ', we have to rewrite the above as follows: We substitute 



r o = r — r + r ■ r" 



(18) 



d 2 W (is) 1 



dudVL (17) 47r 2 c 3 



cti — r+f -r' 



"1 'TCI / [ Q 

jd*r> jdr> lh { , [r x 



x r) 



iuJ(r'^-\-r — r■r , ) 



(19) 



which, if we restrict ourselves to the special case h — 1, means 



d 2 W 



1 



duodVt 4ii 2 c 3 



d 3 r' 



cti— r+f -r' q 

dr' f x -g-jj(r' ,r") x f e ~^ i(r o +r ^ r ' T " ) 

cto— r+f -r' 



(20) 



This is still difficult to evaluate because of the variable bounds of integration. 
We expand them as needed and compensate the error in the integrand: 



d 2 W (19) 1 



duodVt 4tt 2 c 3 



ct' 



Jd 3 r'jdr' g(r' ,r')(r x j A j(r' Q , r' ) x 



R3 ct' Q 



ct'r. D =' min (ctn — r + f ■ f') = const. 

Vr'eK 3 

ct\ D = max (cti — r + r ■ f') = const. 

Vr'eK 3 



(21) 

(22) 
(23) 
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g (r' Q , f ') D = ] jh^^Z fL \ &(r' -cto + r-f- r')G(-r' + cti - r + r ■ f ') 
where Q(x) 




1 > 0, 
otherwise 



Here, the time interval to ■ ■ ■ ti is the detection interval, and t' Q . . . t[ are the 
times for which we have to supply the current data. Thus, the latter mark 
the start and end of our numerical simulation. In the exponent, ^ is large 
and constant, and can therefore be omitted. In fact, it should be omitted for 
numerical reasons: The functions sine and cosine do not work well for large 
arguments. The inversion of (22) and (23) tells us the detection time limits: 



ctn = ct' n + r— min (r ■ f ') (25) 
cti = ct\ + r — max (r ■ f ') (26) 



Let V ro = supp(j(ro,r)) be the support of the current at the time r and 
V — UK the sum of all these. Then, in (21), it is possible to replace the 

Vro 

spatial integral over M 3 by another one over just V, and to calculate the times 
t and ti in (25) and (26) as well (use Mr' G V instead of M 3 ) - at least in 
principle. As we do not want to run our code twice - once to exactly evaluate 
V, the second time to actually calculate the spectrum (to save all currents 
from the first run would definitely exceed any storage capacity) - we simply 
choose V based upon classical estimates of the particle trajectory as a spherical 
volume (which facilitates arbitrary observation directions) of radius R, which 
is large enough to enclose all V ro . This is shown in Figure 1. If the current is 
not localized at the origin, as in the above derivations, but at the middle tk 
of the spherical volume V, then all positions vectors r, r ' have to be moved 
accordingly: 

r^r — rx, r — tk, r= — ^— — — , r=\r\^\r — tk\ (27) 

\r | \r — tk\ 

If pE^y points from the mid point f K towards the observer r, then one obtains 
from (25) and (26), taking into account (27) 



ct = ct' + \r-r K 
cti = ct[ + \r — tk 



- min f — T -^- ■ (f 1 - r K ) ) = ct' Q + If - r K \ + R 
Wev\\r-r K \ ') u 1 

— max — — — r ■ (f ' — r K ) = ct\ + If — f K \ — R 
W'ev\\r-r K \ J 



(28) 
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Changes according to (27) have to be applied to all equations (21) to (24), of 
course. The following explanation is meant to improve the readers intuitive 
understanding of the above purely mathematical description: If we designate 
the edge of volume V that lies furthest away from a given, fixed observer 
by the label a) and the other, closest edge by b), as shown in Figure 1, then 
obviously the time that the radiation needs to travel from a) to the observer is 
longer than from b). The simulation starts at r' = ct' . One would be inclined 
to activate the detector at t = t' Q + [travel time b) to the detector] and to 
close it at t\ = t[ + [travel time from a) to detector] in order to avoid letting 
any radiation go by undetected. However, for consistency, then one would also 
have to take into account any radiation that reaches the detector at time to 
as defined before, but which originates from position a) at a time before the 
beginning of our actual calculation, and about which therefore, no information 
is available. Thus, our conservative approach is to open the detector a little 
bit later. The choice to = t' + [travel time from a) to the detector] throws away 
some radiation of origin b), but now everything that is detected is properly 
taken care of. Similarly, at time t\ as defined above there would be a lack of 
any information about radiation that is coming from b) and superimposing 
the radiation that is still arriving from a). Therefore, the conservative choice 
t\ = t[ + [travel time from b) to the detector] has to be made here, too. 
Considering equation (28), one easily notices that the detection interval can 
become rather short for large values of R. Thus, a spatially extended motion 
can only be analyzed if it is slow, or more exactly speaking, if t' Q and t[ are 
widely separated in time. At least for a particle initially at rest, t' can be 
moved into the past as much as required without altering R, whenever this is 
not the case. 



3.2 Application to quantum mechanics 

In relativistic quantum mechanics, the probability current density is given in 
the following way [i^l 




c-0'7 7^, where 7 = /3 and 7 = (3a, 



(29) 



and can be easily decomposed into the vectorial current density 



Jprob = c^aip 



(30) 



and the scalar density 

Pprob = ^V- 

For the charge current Jin the context of Dirac theory, we choose 



(31) 



3 = Q 3 P rob = qc^ai/j 



(32) 
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This choice is less obvious than it seems. For a extensive discussion of this 
matter, see [jjj ]. The time derivative = ^^Jis derived in the following 

way: 



d _ (32) _ f ^ 

-j = qc—a^ + qc^a- 

9 4 (B = 1} -Aca-i-V- q -l) i> + f3mc^ + qA^ 
at in \ \i c J , 



(33) 
(34) 

(35) 



Note that (34) is an immediate result from the Dirac equation (B.l). The 
momentum operator (35) is transformed by means of Fourier transforms into 
a simple multiplication operator. 



3.3 Reduction to two dimensions 



We have to reduce (21) by means of (B.3) to the two-dimensional case: If 
the observer is located inside the plane, i.e. f ■ r' = r ■ f' ± then equations 
(21), (22), (23) and (24) are transformed into 



d 2 W (2i) 1 



dujdVt 4tt 2 c 3 



ct', 



d 2 r'Jdr' g(r' ,r' ± 



r x 



ct' 



iuj(rQ+r — f-r ^) 



./ (22) . ( , .a -*/ \ 

ct n = min(ct — r + r ■ r , ) 



Vr 



./ (23) ( A v 

cci = max\ct\ — r + r ■ r , ) 

Vf", 



@( — r o + rfi — r + f ■ r ' 



(36) 

(37) 
(38) 

(39) 



To evaluate the above, we need 



fdr'^-^-j ( = g ykrf| + ^ t(5 ^" 
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a 



dt 



it 



+ 



/dtp 
11 at 



(40) 



and as a part of the above 



7 11 Y dt 



\—{ca-\ -V ± - -A + /3mc 2 + ^ U' 

(B.5) [ V « c / J / 



=RJ(0) 



+ (/3mc 2 + gA - g (a ■ 1)) V'} 



/i(J(0). (41) 



Finally, insert equation (41) (apart from the factor hS(0), which has to be 
omitted, cf. appendix B) into (40) and the result into (36). Together with 
(37), (38) and (39) this would give (for the case h — 1) the final formula used 
in our code. Because of its size, we do not print it again. 



3.4 Numerical limitations 



Tests showed that the derivative of the current as shown before in equations 
(40) and (41), although correct in theory, cannot be numerically calculated 
correctly this way. Apparently, the reason is to be found in equation (34): 
This equation is only true, if if) is a exact solution of the Dirac equation. 
Small errors in, for example, the phase of the wave function can have a large 
influence on the result. However, the current (32) itself is independent of the 
phase. Therefore, a less elegant, but numerically superior method is a finite 
difference scheme like the following: 

A higher order scheme would be desirable, however this would require a rather 
impracticable multiple storage of the wave function. Fortunately, At has to be 
chosen only small enough as to let uj max = cover the whole spectral range 
that is to be expected. It can therefore be much larger than the step size of 
the numerical propagation of ip. 
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3.5 Non-relativistic limit 



Just to gain some trust into the above derivation, we now look at its non- 
relativistic and spin-less limit. It can be obtained by ignoring small retardation 
effects r ■ f ' « and using the Schrodinger charge current 4 (which depends 
on the Schrodinger wave function tp t £ C) 



h 



0/m -(^(Wt) - (w t )Vt) - 



(43) 



Similar to the classical theory's limit in subsection 2.2, we insert r ■ f ' « 
and (43) into (20) and, after some work, arrive at: 



d 2 W (20) q 2 



dudVt (43) 47T 2 C 5 



ctlrr (Id 



x r e ■= 



ct -r 



(44) 



Here Tq = ct' and p designates the momentum operator 



pifj t = -ihVtpt Aip t , 

c 



(45) 



whose expectation value is defined, for the case of perfect conservation of 
normalization (iptiipt) — 1 Vt, as follows: 



(46) 



At this point, we can employ the Ehrenfest theorem and, for the case of a 
vanishing magnetic field B « 0, we write 



<9 / Ehrenfest / f^/-> ,\\ 



(47) 



to arrive at a simple formula: 

d 2 W/ (44) q 2 



duodVt (46), (47) 4tt 2 c 5 



cti — r 

dr' Q f x (^—E(x,t')) x f ( 



cto-r 



Vv 



(4* 



Calling a ^ ^.E and (o)±(f) r x (o)^/ x f we obtain 



dr (a)±(f) e" 

eto — r 



(49) 



4 The non-relativistic limit for a particle with spin could be obtained from the Pauli 
theory by insertion of the Pauli charge current, which differs from (43) by an extra 
spin-dependent term 
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27TC 3 



1 

2tt 



ti—r/c 



dt' (a)_ 



-(f) 



-tLUt' 



to-r/c 



to— >— oo 



2ttc 3 



.(50) 



If we expand the integration interval to infinity, as indicated above, and re- 
place a_|_ <-» (a)j_, then we recover the classical result (4) in (50). It returns 
the spectrum as proportional to the Fourier transform of the perpendicular 
component of the mean acceleration. 5 If however, we neglect the retardation 
r in the integration limits of equation (50), then we obtain the result from 
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4 Spectra of free electrons in a laser field 



In order to appreciate the work carried out in section 3, we have made a com- 
parative study of five different methods. In this study, the physical parameters 
are primarily chosen in order to magnify the differences in the spectra and to 
speed up calculations. 

The scenario is the following one: A free electronic wave packet is propagated 
over 20 cycles, two of which are dedicated to the field amplitude's sin 2 -shaped 
turn-on and turn-off phase, respectively. In order to see a reasonable number 
of harmonics in the spectrum of an otherwise free electron 0], such a large 
number of cycles is necessary. A free electron is favourable to the numerics: 
There are, at least in polarization direction, no large canonical momentum 
components, and the wave packet is not torn apart by scattering in addition to 
its unavoidable natural spreading. Therefore, it is possible to use small grids 
of rather low spatial resolution. The large number of cycles implies a high 
frequency in order to keep the overall propagation time short. Note that the 
time resolution for Dirac calculations is almost independent of the problem 7 
and only given by At < -| ~ -Aj. Therefore the frequency was chosen as 
uj = 7.120 a.u. (corresponding to 6.4 nm). 8 In order to stay in the strongly 
relativistic regime - it is only here that we can expect appreciable differences 
- the intensity has to be set to the challenging value E = 300 a.u. (/ = 
3.15 x 10 21 W/cm 2 ). The following calculations were carried out: 

5 Note that most authors ignore the correct prefactor and give the spectrum in 
"arbitrary units" only. 

6 The authors of [3 employ the length gauge E(t) = — VAq with Aq = —E(t) ■ x 
and the dipole approximation A = 0. 

7 Extremely high kinetic energies in the MeV regime are an exception to this rule, 
see |23l ]. They require an even higher temporal resolution since E = -jmc 2 and 7 > 1. 

8 This particular frequency is the highest one that is indicated for the future FEL 
of the so-called „TESLA-Test-Facility (Phase 2)" at DESY in Hamburg Q. 
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the classically-relativistic trajectory of a point particle according to the 
procedure described in section A. 

the quantum mechanical propagation of a Gaussian wave packet according 
to the Dirac equation as shown in . 



the quantum mechanical propagation of a Gaussian wave packet according 
to the Schrodinger equation. 9 



From these data sets, the following spectra were produced: 



• Based on the classical trajectory data, the following was calculated: 

• the relativistically correct spectrum according to equation (1). This will 
be labelled "classical theory (exact)" in the figures to follow. 

• the non- relativistically approximated spectrum according to equation (4), 
labelled "classical theory (approx.)". 

In both cases we used the same (relativistically correct) trajectory data set, 
i.e. the distinction relativistic vs. non-relativistic only refers to the way that 
the spectrum is evaluated. 



• From the single data set according to the Dirac theory, two different kinds 
of spectra were extracted: 

• the fully relativistic one ("Dirac") according to the method presented in 
subsections 3.1 to 3.4: It contains no non-relativistic approximations and 
takes into account all retardation and spin effects. 

• a half-relativistic one ("Dirac-Schrodinger"): the expectation value of the 
electric field was evaluated from the density distribution, and then, ac- 
cording to the procedure in subsection 3.5, the spectrum was calculated. 
In other words, the relativistically correct Dirac probability density was 
fed into a non-relativistic procedure, which was originally developed for 
the Schrodinger probability density. Spin effects are taken care of for the 
propagation, but not for the spectrum calculation, i.e. for the charge cur- 
rent. 

• From Schrodinger data set ("Schrodinger"), the spectrum was evaluated 
according to the method in subsection 3.5. This represents an entirely non- 
relativistic version, with no spin effects at all. 



Additionally, the observation angle $ is varied from 0° to 90° in steps of 
30°. The "Dirac-Schrodinger" method allows us to differentiate between ef- 
fects caused by different wave function propagation methods and those that 
originate from the different treatment of the emitted radiation. 
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4-1 Comparison of the results 



We point out the following features of the spectrum that was calculated using 
the Dirac charge current in Figure 2 and the magnification thereof in Figure 3: 
For $ = 0°, there are only harmonics of first, third and fifth order. Third 
and fifth order are very weak and hard to distinguish from the oscillating 
underground. No red shift is found. At $ = 30°, there are more harmonic 
orders: Up to the seventh order both even and odd harmonics can be seen, 
and a small red shift is found. This effect is enlarged for $ = 60°. For d = 
90° we find the largest red shift and harmonics of up to ninth order. The 
intensity of the fundamental line decreases to a level even below the one of the 
second to fourth order harmonics. The underground intensity also decreases. 
The underlying reasons are best understood by looking at the similar, but 
more intuitive classical picture. Looking at the spectrum in Figure 4 and its 
magnified version Figure 5, which are based upon the classical trajectory and 
the corresponding relativistically correct radiation formula, we notice that for 
$ = 0° there are only odd harmonics, but both even and odd ones up to about 
19th order for all other angles. For order above 20, the peaks disappear in 
the noise. There is no red shift at i? = 0°, but an increasing one for larger 
angles. At = 90°, the intensity of the fundamental line drops below the one 
of the second to fifth order harmonic. The similarities between this classical 
calculation and the Dirac based one are obvious. The main difference is that 
due to the smeared out charge distribution in the quantum mechanical Dirac 
calculation, some sharp structures of the classical spectrum are averaged out, 
especially its higher frequency part. 

The growth of the red shift with increasing angle can be understood qualita- 
tively quite easily (cf. 0): In the rest frame of the electron, the laser frequency 
is red shifted because of the forward drift of the electron in the laser field, and 
the same is true for the oscillation frequency of the electron in polarization 
direction due to the laser's electric field component. If the electron drift is 
directed exactly towards the observer = 0°), the latter sees this actually 
red-shifted frequency in the electron's radiation, which however is Doppler- 
shifted back to exactly its original position in the frequency domain. If the 
observation angle is larger > 0°), the original red shift stays the same, but 
the compensating Doppler blue-shift decreases with increasing angle, leaving 
behind an increasing net red shift. 

The appearance of even and odd harmonics can be understood by the means 
of a symmetry consideration. One needs to consider the approximative for- 
mula (4), i.e. the magnitude squared of the Fourier transform of the acceler- 



9 This calculation was contributed by Andreas Staudt, who is the author of the 
required Schrodinger-code. 
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ation component a± perpendicular to the direction of observation. Reference 
j3] provides the classical analytical solution for the particle velocity j3{rj) in 
terms of the laser phase T]. If we insert a simple harmonic cw-laser field in 
those formulas and take the time derivative, we obtain for the acceleration 
something like 



d - 

a = c—(3(ri) = uir])^ + g(ri)a 2 (51) 

where u(rj) is a rather complicated function, with odd symmetry with respect 
to both rj = und r] = ~, whereas g(rj) features even symmetry with respect 
to r\ = and odd symmetry with respect to 77 = |. In the above, <?i and c?2 
are unit vectors in propagation and polarization direction, respectively. The 
properties of the Fourier transform then automatically lead to even harmonics 
for the <Ti-component of acceleration and odd ones for its ^-component. The 
latter ones are exclusively visible for the observation direction designated by 
$ = 0°, and only the former ones for = 90°. At any angle in between, one 
would expect the appropriate mixture of both even and odd harmonic orders. 
The above analytical discussion explains the complete disappearance of even 
orders in the numerical calculations for = 0° just as well as the suppression 
of the fundamental (which is odd by itself) at i9 = 90°. All deviations from 
the mathematically ideal behaviour have to be attributed to the relativistic 
corrections, the retardation r — f ■ f ' 7^ and last but not least the usage of 
a time-limited, pulse-shaped laser pulse. 

Summing up the above, we can say that the classical and the quantum me- 
chanical description are in good agreement and both deliver well understood 
spectra. One might therefore ask immediately for the motivation of the com- 
putationally expensive Dirac calculations. In fact, this particular scenario was 
set up mainly because we wanted to verify the correctness of the procedure 
described in subsections 3.1 to 3.4 using a system from which one would not 
expect large quantum effects, so any deviations from the well-known classical 
behaviour would easily reveals errors in the method or the code. The agree- 
ment between classical and quantum theory is likely to disappear however 
for laser-driven bound systems. The parameter regime usually referred to as 
"tunnelling regime", for which a description of the "High Harmonic Genera- 
tion" process according to the recollision model [20] is applicable, would be 
of high interest. Unfortunately, this regime is characterized by low frequen- 
cies and therefore long computational runtimes, which renders it impossible 
to study with our Dirac code, at least when simultaneously employing field 
strengths that are large enough to render relativistic and spin-related effects 
visible. Low intensities would be feasible, but not very interesting, as there are 
already a number of completely sufficient Schrodinger 1^1 and Pauli [H, E(j 
calculations for this regime. 
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However, we note some interesting discrepancies when we compare the method 
from subsections 3.1 to 3.4 and the ordinary procedure as shown in subsec- 
tion 3.5. The latter, or more exactly equation (48), can, although originally 
deduced from Schrodinger theory, also be evaluated using the Dirac probability 
density: This is achieved by Fourier transforming the perpendicular acceler- 
ation component a± as in equation (50), where the acceleration is taken as 
the expectation value of —E(x, t) with respect to the Dirac wave function. 
As an example, this is shown in Figure 6 for an observation angle $ = 30° 
and enlarged in Figure 7 (second graph from top to bottom in both figures) 
side by side with the exact Dirac procedure, the also exact classical method 
according to equation (1) and its non-relativistic approximation (4), as well as 
the non-relativistic Schrodinger calculation 10 , which is also based upon equa- 
tion (50), but is derived from a non-relativistic Schrodinger wave function. As 
one can see, the three not entirely relativistic methods ("Dirac-Schrodinger", 
"classical theory (approx.)", and "Schrodinger") significantly over-estimate 
the red shift 11 , because in these calculations, the Doppler effect due to the 
electron's drift in propagation direction is not included. The even harmonics, 
which should be visible for any non-zero observation angle, are missing in both 
spectra according to equation (50) ("Dirac-Schrodinger" and "Schrodinger"). 
This is because only the perpendicular component of the electric field is taken 
into account, weighted either with the Dirac or the Schrodinger probability 
density. As a result, the spectrum features a simple cos($)-dependence and 
even vanishes for $ — > 90°. There is no change in structure with varying angle, 
it has always the shape of the version for t? = 0°, although it is only valid for 
this particular angle. 

Apart from the methodical differences laid out above, we also see differences 
between the two classical calculations and all three quantum mechanical ones. 
Especially the higher order harmonic peaks that are visible in the classical 
curves are averaged out in the case of the quantum spectra. 



Note: The Schrodinger calculation only seems to exhibit far less noise than all 
the other calculations. The reason is that it was carried out using a much coarser 
frequency resolution (about one fifth). Apart from that, its prefactor is not known. 
To indicate this, the vertical scale in the corresponding plots was omitted. 
11 Note: These three methods even predict a red shift at d = 0°, which is not shown 
here. 
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5 Conclusions 



From the above, it can be concluded that the procedure laid out in detail 
in subsections 3.1 to 3.4 is superior to the established techniques in terms of 
angle dependence of especially the even harmonic orders, and in terms of the 
predicted red shifts, whenever the parameters are chosen in a way that leads 
to a significant influence of retardation and magnetic field effects. In principal, 
even spin effects 12 could be analyzed this way, since the method is based on 
the Dirac theory. Although originally developed for the usage in combination 
with Dirac wave functions, the method could be employed in the context of 
the Schrodinger or Pauli theory as well 13 , if one abandoned the Ehrenfest 
theorem and evaluated equation (44) directly, and it would at least correct 
the angle dependence. Employing equation (20) would even lead to correct 
red shifts, and the only remaining quantitative deviations would be due to 
the underlying non-relativistic wave packet propagation. The computational 
cost for the spectrum alone is not essentially higher than for the more sim- 
ple schemes, therefore its general usage can be recommended. However, the 
parameter regime in which its advantages are best visible are already very 
demanding as far as the wave packet propagation is concerned. 

For the free electron scenario, classical results are very similar to the ones ob- 
tained with our new method. We see no true quantum effects apart from the 
disappearance of the higher order harmonics. This may be attributed to aver- 
aging effects by the spatially extended Dirac charge current. More significant 
effects are very likely to be observed once that atomic systems are considered. 
It must be noted, though, that our main goal was to devise a relativistically 
correct scheme, as opposed to the non-relativistic existing ones. Our results 
show that the relativistic corrections are significant for the parameters em- 
ployed here. 
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We have chosen an initial state that has spin up with respect to the magnetic 
field direction 173, but did not analyze the influence of variations of the initial spin 
polarization on the spectrum. Spin effects such as those studied in j^H] would even 
require three-dimensional calculations. 

13 Note that in the Pauli theory, the charge current would have to be corrected for, 
see subsection 3.5. 
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A Classical equations of motion 



In classical relativistic mechanics, momentum change equals the applied force, 

(A.l) 




where the force / is the well-known Lorentz-force: 



fej t t) = qE(gj) + ^vx B^ >t ). (A. 2) 



Equation (A.l), which represents a second-order differential equation with 
respect to the three-dimensional position variable x can be easily transformed 
[l0( into a first-order system of differential equations in terms of the six- 
dimensional combined variable (x,v): 



j t *= y i - ^ (fmt) - i (V • f W )) v) (A.3) 

^-x = v (A.4) 
dt v ' 

Formally, the above means 

^-(x,v) = G^^t) ( A - 5 ) 

where Gig^t) represents the six- dimensional function that is obtained from the 
right hand side of equations (A.3) and (A.4) after inserting the force (A. 2). 
The first-order system (A. 5) is then easily numerically integrated by means of 
a standard fourth-order Runge-Kutta algorithm j36|. 



B The Dirac equation 



The Dirac equation in 3 + 1 dimensions reads as follows: 

ih^ = col ■ \ -V - - A ) ib + (3mc 2 ip + qA ib. (B.l) 
at \% c ) 

In the above, h is Planck's constant, q and m represent the particle's charge 
and rest mass, c is the speed of light, a 1 and j3 are the usual Dirac matrices 
(i = 1, 2, 3) [22|. We use atomic units, so for an electron h — m — — q — 1 a.u., 
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and c = 137.036 a.u. Suppose we have a potential that is invariant under 
translations in one particular direction that is designated by n 



A = A (x±), A = A(x±), where x± = n x (x x ft), (B.2) 
then we can use the separation ansatz 

vf)(x) = i/j'(x±)x(x\\) an d x( x \\) — ex P (^l^ll) ' wnere x\\=n-x 

(B-3) 

in the Dirac equation (B.l). The result is: 

dip' ( h -> q -\ 
ih-Qj- — ca-i-Vx +P\\n A j ip' +j3mc 2 ?p' +qA ip' , where V_l = nx(Vxn). 

(B.4) 

In general, there will not be a single momentum pjj, but rather a momentum 
distribution giving specific weights for arbitrary p» e M and the solution is a 
wave packet instead of (B.3). We restrict ourselves to the simple choice pjj = 
to eliminate this constant. Other values would not provide more physical in- 
sight, and true wave packet calculations would cause the computing times to 
grow into the regime of true three-dimensional calculations. Note however, 
that the ansatz (B.3) inevitably leads to a meaningless factor M(0) in bilinear 
expressions of such as expectation values or the integrated current density. 

We finally obtain an effective Dirac equation in two dimensions: 

= C 5 . [ ^v_l - -A] if/ + (3mc 2 i)' + qA ij' (B.5) 
ot \l c J 

The above equation is integrated 01 using the split-operator scheme [38^ 
presented in our earlier work . The dynamical grid techniques described in 
there are especially well-suited for the free electron problem discussed here, 
and dramatically reduce computation times. 

A vector potential that describes a laser beam of electric field amplitude Eq 
and frequency lu, which is linearly polarized in direction a and propagates in 
direction of the wave vector k, such as 



A= -a—E gfoifut- k ■ x), A = 0, (B.6) 

UJ 

is invariant in direction n T> =' k x a. As can be easily seen, in this context "two- 
dimensional" does not mean that there is no magnetic field B (which points in 
the xn-direction). Also note that the above two-dimensional treatment is exact 
for the special case of a free electron in a laser field that we considered in this 
paper, because the vector potential (B.6) features the required symmetry. 
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Fig. 1. Sketch of the radiation scenario: The three black-shaded areas that carry time 
(t') labels represent the electron at the beginning (t' = t' ) and the end (t' = t^) 
of the simulation, as well as at an arbitrary point in time (t' 6 [io>£i]) in between, 
each enclosed inside the dashed boundaries of the co-moving and -growing numerical 
grid. The oscillatory curve is the trajectory of the particle's centre. The observer 
(or detector), which is activated during the time interval to • • • t±, is located under 
an angle with respect to the laser propagation direction at a distance \r — fk\ as 
measured from the middle of the volume V (large grey-shaded circle), which is a 
sphere of radius R and origin tk- At the same time, the radius R represents half the 
diagonal line of the solid rectangle, which is the enclosing envelope of all occurring 
numerical grids. The two parallel lines perpendicular to the line of observation 
labelled "a)" und "b)" mark two specific distances to the observer mentioned in the 
text. 
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Fig. 2. Spectra according to the Dirac calculation: The plot shows the spectrum, 
evaluated according to the method from subsections 3.1 to 3.4, of a free electron 
that interacts with a laser pulse of 20 cycles (two of which are reserved for turn-on 
and turn-off) with an amplitude Eq = 300 a.u. (/ = 3.15 x 10 21 W/cm 2 ) and fre- 
quency u) = 7.12 a.u. (A = 6.4 nm). The fully relativistic spectrum due to the 
two-dimensional Dirac charge current is shown for several different angles of obser- 
vation $ with respect to the propagation direction. 
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Fig. 3. Spectra according to the Dirac calculation; zoom of Figure 2: The appearance 
of the even orders for non-zero angles and the increase of both intensity and red 
shift of all harmonics with growing observation angle is clearly visible. At $ = 90°, 
the fundamental line is suppressed. 
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Fig. 4. Spectrum according to the classical calculation: The plot shows the spec- 
trum of a free electron in interaction with a 20 cycle laser pulse (including two 
cycles for turn-on and turn-off, respectively) with an amplitude of Eq = 300 a.u. 
(I = 3.15 x 10 21 W/cm 2 ) and a frequency of u = 7.12 a.u. (A = 6.4 nm). The fully 
relativistically calculated spectrum due to two-dimensional motion of a classical 
point particle is shown for several different angles of observation -d with respect to 
the propagation direction. 
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Fig. 5. Spectrum according to the classical calculation; zoom of Figure 4: The ap- 
pearance of the even harmonic orders for non-zero angles and the increasing red 
shift of all harmonics with growing observation angle are easily visible. At •& = 90°, 
the fundamental line is suppressed. 
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Fig. 6. Comparison of several methods to evaluate the radiation spectrum: The 
top graph is based on the method from subsections 3.1 to 3.4 and shows the fully 
relativistic spectrum due to a two-dimensional Dirac charge current. The second 
graph from top to bottom displays the spectrum obtained by Fourier transforming 
the perpendicular acceleration component a±, where the acceleration is taken as 
the expectation value of ^-E(x, t) with respect to the Dirac wave function, similar 
to the way shown in subsection 3.5 for a Schrodinger wave function. The third 
graph shows the fully relativistic spectrum calculated for the classical trajectory 
according to subsection 2.1, and the fourth graph the one that can be obtained from 
the same trajectory data in the non-relativistic limit according to subsection 2.2. 
The bottom graph displays the result according to the method in subsection 3.5 
for a non-relativistic Schrodinger calculation. No vertical scale is given for this last 
calculation, which is a contribution of Andreas Staudt, because the prefactor is not 
known. 
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Fig. 7. Enlarged view of Figure 6: The differences of the distinct methods in terms of 
red shift and harmonic orders are obvious. All methods that contain non-relativistic 
elements (Dirac-Schrodinger, approx. classical theory and Schrodinger) feature sig- 
nificant red shifts at the observation angle displayed here (30°). For the exact classi- 
cal and the Dirac-calculation noticeable red shifts only occur for much larger angles. 
Apart from that, even harmonic orders are missing in both the Dirac-Schrodinger 
and the Schrodinger calculation. The latter only seems to have less numerical noise 
than all the others, but this is due to a coarser frequency resolution. 
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